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ABSTRACT 

Context. The first million years of the pre-main-sequence (PMS) evolution of low-mass stars are characterized by magnetospheric 
accretion, a circumstellar proto-planetary disk, and the processes leading to its dispersal. Among these, photo-evaporation caused by 
strong X-ray emission from the central star is probably significant. Several aspects of the X-ray emission from coronae and accretion 
shocks remain mysterious, e.g., whether and how much accretion affects coronal emission. 

Aims. We studied the X-ray variability of ~lMyr old low-mass PMS stars as a function of timescale, stellar rotation, and stellar 
characteristics to gain insights into the working of PMS coronae, their X-ray emission, and the circumstellar environment in which 
they are immersed. 

Methods. We have exploited the ~850ksec long Chandra observation of the Orion Nebula Cluster obtained by the COUP collabora- 
tion in January 2003, and statistically analyzed the X-ray lightcurves of low-mass stars in several subsamples. Our main focus was to 
characterize the different X-ray behavior of stars with and without circumstellar accretion disks, and to infer the physical mechanism 
responsible for the observed variability. 

Results. Accreting stars (classical T Tauri stars, CTTSs) are found to be more variable than non-accreting stars (weak-line T Tauri 
stars, WTTSs) at all timescales and in all X-ray energy bands considered. Variability is seen to increase with timescale, up to the 
longest probed, ~10 days. Signatures of rotational modulation are observed for both CTTSs and WTTSs, and most clearly for CTTSs 
in the soft X-ray band. Lower mass stars are more variable than higher mass stars. 

Conclusions. We propose that the difference in variability between CTTSs and WTTSs may be explained assuming that the X-ray 
emission of CTTS is affected by time- variable absorption due to circumstellar structures, such as warps in the inner disk and/or accre- 
tion streams. This suggestion is appealing because, in the hypothesis that the coronae of CTTSs and WTTSs are similar, it may also 
explain why CTTSs have lower and more scattered X-ray emission levels compared to WTTSs. 

Key words. Stars: activity - Stars: coronae - Stars: variables: T Tauri, Herbig Ae/Be - Stars: formation - Accretion, accretion disks 
- X-rays: stars 



1. Introduction 

The evolution of young stars in the first few million years after 
they emerge from their protostellar dusty envelopes is character- 
ized by the presence of circumstellar disks and intense magnetic 
fields. These latter crucially regulate the interaction between the 
star and the disk, they mediate for instance the angular momen- 
tum transfer between the two and channel streams of accret- 
ing material from the disk onto the stellar surface (e.g. |Bouvier| 
|et al. 2007b[ l. This material is shocked in the impact, produc- 
ing optical, UV, and soft X-ray "excesses" ([Calvet & Gullbriiig 



1998 ). Magnetic fields are also responsible for the confinement 



and heating of X-ray-bright coronae (Feigelson & Montmerle 
1999[ l, whose high energy radiation, along with that from accre- 



tion shocks, is likely to significantly influence the evolution of 
circumstellar disks: soft X-rays ionize and heat parts of the disk 
and may be the main cause of disk photo-evaporation and, con- 
sequently, of their eventual dispersal (e.g. |Glassgold et al.|2004[ 
?). 

The above processes are highly dynamical. Strong variabil- 
ity at all wavelengths and on many different timescales is indeed 
a distinctive trait of young accreting pre-main-sequence (PMS) 
stars Poy|1945| |Alencar et al.|2010| , also known as classical T 



Tauri stars (CTTS). Weak-line T Tauri stars (WTTS), similarly 
young stars for which accretion signatures are not observed, rep- 
resent almost certainly the following stage of PMS evolution, 
once the circumstellar accretion has stopped and the dusty disk 



1981 



has dissipated or is quickly dissipating (jFeigelson & Kriss 
Lada 19 87). Compared to CTTSs, WTTSs are much much less 
variable in the optical band and their lightcurves show mostly 
periodic variations attributable to the modulation of cool stellar 
spots ( Alencar et al.|2010 1. 



The mechanism responsible for the strong and irregular opti- 
cal variability of CTTSs, although likely related to accretion, re- 
mains elusive. The two main hypotheses are /) variability of the 
emission from the accretion shock(s), both intrinsic and caused 
by rotational modulation, or //) variable absorption of the photo- 
sphere, or part of the photosphere, from unstable and optically 
thick accretion streams and/or warps in the circumstellar disk. 
Grankin et al. ( 2007| 2008| l, for example, suggested from a sta- 
tistical study of CTTS that in about 25% of the cases variable 
absorption is the likely cause of optical variability, while time- 
variable accretion is favored in the remaining cases. However, 
owing to the similar effects of spots and absorption on the broad- 
band lightcurves, the absorption scenario could not be excluded 
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in most cases. The thorough study of the CTTS AATau by 
Bouvier and collab orators (Bouvier et al."1999', '2^3] |Menard| 
let al.|2003| |Bouvier et al. 2007a; Grosso et al. 2007| Tias shown 
that for this near edge-on star-disk system, the strong, irregu- 
lar optical variability is readily explained by occultation of a 
significant fraction of the stellar surface by a warp in the in- 
ner disk. This warp, located at the corotation radius, produces 
quasi-periodic eclipses with the period of the stellar rotation. 



Alencar et al.|(2010l, using high-quality optical lightcurves of 



stars in the NGC2264 region obtained with the CoRoT satel- 
lite, have recently proposed that an AA Tau-like scenario might 
indeed explain the optical variability of at least ~30% of the 
CTTS. The A A Tau-like phenomenology was also confirmed by 
IMorales-Calderon et al.| ( |20n] l, from Spitzer/IRAC (3.6/4.5//m) 
lightcurves of ONC sources, with a much lower incidence, ~5%, 
than found by |Alencar et al. ( 2010 1, however. 

In the X-ray ba nd all PMS stars, both CTTS and WTTS, 
are intense sources ( Feigelson & Krissl|1981 Flaccomio et al 



|200 3b ; Preibi sch et al.||2005| l. Thermal emission from magnet- 
ically confined plasma (i.e. coronae) at ~ 10-30 MK accounts 
for most of the observed emission. Accreting stars (CTTS) also 
show a softer plasma component at ~2-4MK that, given its 
temperature and plasma density, can be interpreted as the di- 
rect product of the accretion shock. This component has only 
been unambiguously detected in the few CTTS for which high- 
resolution X-ray spectra were obtained (e.g. Kastner et al._2002, 



Giidel & Telleschi|2007 Argiroffi et al.|2011[ l, as demonstrated 
by emission-lines from cool plasma and by density-sensitive He- 
like triplets. It is natural to think that such a component is a com- 
mon feature of CTTS. However, for most stars, with the possible 
exception of TW Hya ( |Kastner et al.|2002|[Dupree et al.|20T2l l, 
this cool plasma does not contribute significantly to the X-ray 
flux detected in the 0.3-lOkeV energy band covered by present- 
day imaging detectors, such as those onboard of Chandra and 
XMM-Newton. In fact, in spite of this additional soft component, 
CTTS are, on average, two to three times fainter in X-rays com- 
pared to WTTSs of the same mass or Lboi ([Neuhaeuser et al 



1995 



Damiani et al.|[T995| [Flaccomio et al.||2003a[ [Preibisch 



et al. 2005). Moreover, their activity levels at any given mass 
or Lboi show a larger scatter (Preibis ch et al.[ |2005 ). These dif- 
ferences cannot be ascribed to differences in stellar rotation, the 
main factor that determines activity levels on MS stars, and re- 
main largely unexplained. A higher average absorption of CTTS 
with respect to WTTS, which would lower the observed X-ray 
fluxes of CTTS (but not the intrinsic ones), was excluded by 
[Preibisch et al.[ ( |2005| l on the basis of the Chandra Orion Ultra- 
deep Project (COUP) data (Getm an et al. |2005 ), which allowed 
determining individual absorption corrections for all stars. These 
authors, as well as others, instead proposed that accretion might 
modify the magnetic field structure (and thus coronae, see also 
[Flaccomio et al.|2003b[ [Telleschi et al.|2007 ^ or the stellar struc- 
ture in such a way as to influence the stellar dynamo. [Telleschi] 
et al. ( 2007[ l in particular, based on data from the XMM-Newton 
Extended Survey of the Taurus molecular cloud (XEST), pro- 
posed that the accreting material from the disk cools down part 
of the coronal plasma to temperatures outside the X-ray regime. 



Gregory et al. (2007) modeled coronal and accretion structures 



by extrapolating magnetic field surface maps (from Zeeman- 
Doppler imaging) and concluded that the lower observed X-ray 
luminosity of CTTS might indeed be due to the obscuration of 
large parts of the corona by dense gasseous accretion streams. 
This scenario is closely related to that invoked to explain the 
AA Tau-like optical variability, the main difference being that 
for X-rays the absorbing material would be the gas in the ac- 



cretion streams rather than the dust in the inner disk warps at the 
bases of the same accretion steams. The accretion streams would 
only intercept a fraction of the emitting coronal plasma and the 
emerging X-ray flux would thus be composed of two compo- 
nents, only one of which is heavily absorbed. As also discussed 
in this paper, correctly correcting the observed X-ray flux for 
extinction would be almost impossible, at least with the spectral 
resolution and sensitivity of present-day instruments. 

In addition to the differences in observed X-ray flux, other 
differences in the X-ray emission of CTTSs and WTTSs have 
been reported. Flaccomio et al. ( |2006| l, for example, reported 
tentative evidence that CTTSs in the NGC 2264 star-forming re- 
gions are more X-ray variable than WTTS. The same data, as 
weU as the XEST data ( [Telleschi et al.|2007[ ), also indicate that 
the plasma on CTTSs is, on average, hotter than that on WTTS. 
At the same time, some CTTSs show evidence for a plasma com- 
ponent cooler than seen on WTTSs. While this latter observa- 
tion may be attributed to plasma heated in accretion shocks, the 
harder and more time- variable emission could be consistent with 
time-variable absorption from circumstellar material (accretion 
stream and/or disk warps). 

From this discussion it is clear that studies of X-ray vari- 
ability may provide useful information on the mechanism that 
regulates either the observed or the intrinsic X-ray activity level 
of PMS stars. All T Tauri stars are strongly variable in the X- 
ray band on timescales of hours and days: the most prominent 
form of variability is undoubtedly represented by flares, with 
short rise phases during which the flux reaches up to ~ 100 times 
the pre-flare level, followed by slower decay phases lasting sev- 
eral hours (see e.g. Favata et al. 2005J for a detailed study of 



bright flares in the COUP observation; [Wolk et"ar 



2005j and 



[Caramazza et al.[[2007[ for statistical studies). Modulation due 
to the inhomogeneity of coronal structures and stellar rotation 
has also been detected by [Flac comio et aT] ( [2010| l in the ~13- 
days long COUP observation (including gaps; the actual expo- 
sure time was ~850ks). This work also highlighted the difficulty 
in separating flaring from other forms of variability. However, 
variability on timescales longer than one day has generally not 
received much attention, partly because of the scarcity of appro- 
priate datasets. Moreover, the highly stochastic nature of the X- 
ray emission has hindered statistical studies trying to correlate 
the "average" properties of variability with stellar and circum- 
stellar characteristics. 

An interesting result has recently been obtained by us from 
the correlation of X-ray and optical emission. In a previous 
study, [Stassun et al.[ ( [2007 1 had already searched in vain for a 
correlation between the COUP X-ray lightcurves and simulta- 



neous ground-based optical multi-band lightcurves. Flaccomio 
jet al.[ ( |20T0] l have instead employed simultaneous space-based 
observations of NGC 2264 members by Chandra (X-rays) and 
CoRoT (optical, broad-band). By comparing X-ray and optical 
fluxes obtained in two 30 ks time windows separated by 16 days, 
we have found that, for CTTSs only, the X-ray variability in the 
soft 0.5-1.5 keV band is significantly correlated with the opti- 
cal variability. This correlation is seen neither for WTTS nor 
for the harder X-ray band (1.5-8.0keV). We interpret these find- 
ings as caused by the effect of time-variable absorption from cir- 
cumstellar material. From the relative magnitude of the X-ray 
and optical variability we conclude, moreover, that the absorb- 
ing material must be dust-depleted, thus pointing toward accre- 
tion streams within the dust-sublimation radius, as proposed by 
Gregory et aL| ( [2007| . 



If the X-ray emission of CTTSs is indeed absorbed by cir- 
cumstellar structures in co-rotation with the star, their X-ray 
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variability should /) be more extreme than on WTTS, //) reach 
maximum amplitudes on the timescale of the stellar rotation, 
and ///) show clearer signs of rotational modulation. Most X-ray 
observation of nearby young stars performed with present-day 
instruments are shorter than ~100ksec (~ld) and therefore un- 
suitable for studying variability on the timescale of typical stel- 
lar rotations (-1-15 d for stars in the ONC; 'Herbst et al.'2002'). 
The most notable exception is the almost uninterrupted 850 ks 
long COUP observation of the ONC. We have thus decided to 
exploit this dataset once more to investigate X-ray variability as 
a function of timescale from hours to ~12 days, thus covering 
the stellar rotation of most ONC stars. In addition to its length, 
the COUP observation with its ~1600 X-ray sources also pro- 
vides the advantage of a rich and optically well-characterized 
sample of PMS stars. Given the highly stochastic nature of vari- 
ability, large stellar samples are crucial to reliably determine "av- 
erage" variability properties for different stellar subsamples. We 
can then try to correlate these average properties with stellar and 
circumstellar characteristics to investigate the physical origin of 
variability. 

This paper is organized as follows: §[2] describes the COUP 
dataset; §|3]illustrates the two methods we used to quantify vari- 
ability, one non-parametric (§ |3.1| l and one tailored to impulsive 
flares (§ |3.2[ );§[4| describes the selection of the stellar samples 
used for the analysis. Section[5]presents the result of the variabil- 
ity analyses, which are then more fully discussed in §|6] along 
with their implications. Finally, we summarize our main conclu- 
sions in §|7] 

2. The COUP data 



The COUP observation is described in detail by Getman et al. 
( 2005 | l along with the procedures followed for data preparation, 
source detection, identifications with optical/NIR catalogs, and 
photon extraction. We used the products of this analysis and, in 
particular, the original lists of source and background photons 
extracted for each source. We also adopted from Getman et al. 
( 2005| l the stellar and circumstellar parameters for the COUP 
sources, both collected from the literature or derived by these au- 
thors (e.g. magnitudes, masses, etc.). We updated this database 
with respect to stellar rotational periods and classifications in 
terms of near-IR excesses indicative of disks and envelopes. We 
used the same collection of rotational periods as Flaccomio et al. 
( 2005| l, with the addition of 48 new periods derived from an 



-250 h Spitzer observation of the ONC performed as part of the 
"Young St ellar Object VARiability" (YSOV AR) program and re- 
ported by [Morales-Calderon et al. (2011 1. From this work we 



also adopted the classification in terms of spectral energy distri- 
bution (i.e. Class I, Class II, and Class III). 

3. Analysis methods 

We first discuss the non-parametric variability analysis that is at 
the basis of our following discussion, and then the method we 
used to asses the contribution of X-ray flares to the total source 
variability. 

3. 1 . Variability amplitudes 

To study how the X-ray fluxes vary as a function of timescale, we 
subdivided the 850 ks COUP observation into shorter segments: 
we chose to adopt segment lengths of 30 ks, as a compromise be- 
tween being able to probe short timescales, specifically with re- 
spect to the typical stellar rotation periods, and the need to have 




6«10 

Time [sec] 



Fig. 1. Example X-ray light-curve for a bright COUP source. 
The horizontal segments, superimposed on the lightcurve (black 
line), indicate time intervals and relative photon fluxes for the 25 
30 ksec segments on which the variability analysis is based. 



a sufficient number of detected photons per segment and thus 
limit Poisson uncertainties. With this choice we selected 25 time 
segments within the COUP observation, for each of which we 
measured the background-subtracted photon flux. Figure[T]illus- 
trates, superimposed on a binned light-curve (£=0.5-8.0 keV) of 
one of the brighter COUP sources, the time segments we adopted 
for the following analysis. 

In general terms, our analysis consists of quantifying the 
variability between segments separated by a given time inter- 
val. At, by taking the ratio between the fluxes in the first and 
second segment, FijFi = F(t)/F{t + Atj^ We then built dis- 
tributions of these variability measures (i.e. flux ratios) for all 
segment pairs within the same lightcurve separated by At and 
for all X-ray sources in given subsamples. We then study how 
the width of the distributions of variability measures obtained 
for the same subsample varies as a function of timescale. At, and 
compare the trends for different subsample^ For our purposes 
the width of the distribution is defined as the difference between 
the 75% and the 25% quantiles. We refer to this width as the 
variability amplitude iVA) of the sample at the given Af. For 
stars for which we also know the rotation period, we also explore 
the dependence of VAs on rotational phase differences. To this 
purpose we construct and compare the widths of distributions of 
flux ratios taken between segments separated by rotational phase 
differences: = AtjProt- 

The whole analysis was carried out considering three energy 
bands: the full Chandra band (0.5-8.0 keV), a soft band (0.5- 
1.5 keV), and a hard band (1.5-8.0keV). As an example, we plot 
with a thick broken line in Fig.|2] the differential distributions 
of \og{Fxl F2) obtained for the full X-ray spectral range and for 
Af ~8 days (effectively, 7.1< Af <9.6 days). In this example 
the stellar sample comprised 167 Class III stars (see below for 
sample definitions). Given the adopted Af range, each of these 
stars contributes 47 flux ratios and the distribution is defined by 
7849 datapoints. 

The finite width of the log(Fi/F2) distributions, i.e. their de- 
parture from log{Fi/F2) - 0.0, can be caused by two factors: 
real flux variability, or Poisson uncertainties in the number of 
detected photons. To isolate this last unwanted contribution, we 
have taken two approaches: first, we limited our search sample to 
X-ray sources with "good" statistics; second, we quantified the 
effect of Poisson noise on the variability distributions through 



' We experimented with different ways to measure variability be- 
tween segments, all yielding consistent results. 

- The number of datapoints (flux ratios) that define the distributions 
is obviously dependent on timescale, being generally lower at longer 
At. 
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Fig. 2. Normalized histogram of all flux ratios, log(Fi /F2), taken 
between 30ksec time intervals separated by At =7.1-9.6 days 
(logarithmic mean: 8.24 days), for the 167 Class III stars in our 
flux-limited sample (broken line). The narrower distribution, ac- 
tually a bundle of ten almost indistinguishable ones, represents 
the contribution of Poisson statistics to the observed distribution 
as determined from Monte Carlo simulations (see text). 



Monte Carlo simulations. More specifically, we imposed a limit 
of 750 counts within the 25 chosen time segments, i.e. we re- 
quired >30 counts per segment (or 1 count/ksec) on average. 
This condition reduces the number of available X-ray sources 
from the initial 1616 to 496, 319, and 361 for full, soft, and 



hard band, respectively. As discussed in more detail in §6.3 



this selection introduces a bias on the resulting samples, quite 
directly in terms of Lx, and indirectly, through the correlation 
of Lx with fundamental stellar properties (e.g. Preibisch et al. 
2005 ] !, in terms of stellar mass and bolometric luminosity. 



To evaluate the effect of Poisson uncertainties on the flux- 
ratio distributions, and ultimately on the VAs, we generated for 
each of the considered distributions (i.e. for each sample and 
timescale) 1000 simulated flux-ratio distributions assuming that 
sources are intrinsically constant and that the measured fluxes 
only vary because of Poisson fluctuationsj^Figure|2]shows ten 
of the simulated constant-source distributions (the bundle of nar- 
row distributions centered at log(Fi/F2) = 0.0). In this case it 
is quite clear that Poisson uncertainties cannot account for the 
widths of the actual flux-ratio distributions. In general, we quan- 
tified the contribution of Poisson statistics on the VA of each of 
the considered samples and timescales by recording the medians 
and Icr dispersions of the VAs measured from the 1000 simu- 
lated constant-source fluxes. 



3.2. Flare analysis 

Coronal flares, characterized by short rise-phases and slower de- 
cays, are an obvious source of variability for our sources. To 
discuss the contribution of flares to the overall variability ampli- 



^ More specifically, for each pair of time segments involved in the 
flux-ratio distribution we generated simulated source (and background) 
fluxes, from Poisson distribution appropriate for a flux equal to the mean 
of the two measured fluxes. 



tudes (§6.1 1 we seek to quantify flaring activity statistically for 
a given stellar sample. X-ray flares in the COUP datasets have 
already been investigated by at least four studies (Favata et al. 
2005l |Wolk"et al. 2005; Car amazza et al.|2 007; Johnston eetaL 

|W( 



2012 



. We adopted the technique used by [Wolket al., ( ,2005] rand 



Caramazza et al.| ( |2007| l, that is, we detected individual flares 
using an automated procedure and then derived their rate of oc- 
currence as a function of their intensities as well as the distribu- 
tion of "flare durations". A full description of the method can be 
found in [Caramazza et al.| ( |2007| l. 

4. Sample definitions and selection criteria 

As described in the previous section, our parent samples for 
the three X-ray spectral bands comprise all COUP sources with 
an average count-rate >1 count/ksec in the 25x30 ks selected 
time intervals. We then selected physically relevant and homo- 
geneous stellar samples for our study of variability as a function 
of timescale and stellar characteristics. 

First, we excluded high-mass stars, M > 3M0, the X-ray 
emission of which is of non-coronal origin, or comes from a dif- 
ferent kind of corona, or originates in unresolved lower mass 
companions (e.g. jDamiani et al."1994; Zinnecker & Preibisch] 
1994:,Hamaguchi et al.|2005t,Stelzer et al._2009 ). Then, to fur- 
ther restrict the range of masses, and that of absorptions, we ex- 
cluded stars with Ic >16, corresponding to M > 0.1 - O.2M0 for 
COUP sources with low optical absoiptiorj^(Av < 2.0). Finally 
we excluded proper-motion non-members according to the col- 
lection of Hillenbrand| ( [T997| , unless they had indications of cir- 
cumstellar accretion or circumstellar material (see below). With 
these selection our search samples are reduced to 369, 289, and 
261 stars for the full, soft, and hard bands, respectively. 

Next, we defined the following subsamples of low-mass stars 
based on their accretion properties, presence of circumstellar 
disks, and mass: 

- CTTS and WTTS, based on the equivalent width ( EW) of the 
8542 A Ca II line, following the same approach as Flaccomio 
et al.| (|2003b') and ' Preibisch et al.| (|2005|. Namely, we se 



lected as accreting stars (CTTS) those with the Ca II line in 
strong emission, EW<-1, and as non-accreting stars (WTTS) 
those whose Ca II line clearly is in absoption: EW>1. By 
excluding intermediate EW values, this approach results in 
better defined and extreme samples: the selected CTTS, in 
particular, will be biased toward high accretion rates. 

- Class II and Class III, based on the observation of mIR ex- 
cesses due to warm circumstellar material, as reported by 
[Morales-Calderon et"aLl ( |20n] i. 

- Low- and intermediate-mass Class II and Class III stars, de- 
fined as complementing the above selections with the follow- 
ing constraints on mass: M < 0.5Mq and M > 0.8Mq. 

We considered other possible subsamples, but, since their 
sizes are too small to derive meaningful results, we will not dis- 
cuss them. For the samples and X-ray energy bands discussed in 
the following, Table[T] summarizes the size of the sample and, if 
relevant, the size of the subsample of stars for which we know 
the rotational period. Most of our analysis will be based on the 



For the low-mass limit we chose to rely on the magnitude instead 
of the masses estimated from the HR digram and evolutionary models. 
This is because this latter determination of masses, ultimately based on 
spectral types, is less complete at the low-mass end. We would thus 
have lost a number of otherwise useful stars for the following analysis. 
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Fig. 3. Distributions of count rates for our flux-limited samples 
of CTTS (red), WTTS (green). Class II (orange), and Class III 
(blue) stars. Sample sizes and median count rates are reported in 
the legend. Count rates are defined here as the total number of 
counts detected in the 25x30 ksec time intervals under consider- 
ation, divided by the exposure time, 750 ksec. 

Table 1. Sizes of the analyzed samples. Figures in parentheses 
refer to subsamples with known rotational periods 



, All masses <0.5Mo >0.8Mo 
Sample 

^ 0.5-8kev 0.5-1. 5kev 1.5-8kev 0.5-8kev 0.5-8kev 

CTTS 57 (19) 33 (11) 41 (13) - - 

WTTS 110(72) 93 (62) 82 (53) 

Class II 145 (70) 109 (54) 109 (52) 55 27 

Class III 164(105) 138 (94) 113 (73) 73 43 



samples of CTTSs, WTTSs, Class II, and Class III stars, which 
we therefore refer to as our main samples]^ 

In the following we used of the statistical properties of the X- 
ray sources in our samples to constrain the origin of their X-ray 
variability. Figure[3] shows the distribution of (full-band) count- 
rates for our main samples. In spite of our choice of count-rate 
limited samples, and consistent with the results of Flaccomio 
( |2003| ) and Preibisch et al.,(2005) , CTTS have statistically lower 
count-rates with respect to WTTS (a factor of ~2 in the median). 
A difference is also seen between Class II and Class III, but it 
is smaller, possibly owing to the less extreme nature of the se- 
lection criteria. Significant diff'erences in median count-rates are 
also observed between low- and intermediate-mass Class II and 
Class in stars, with the higher mass samples being systemati- 
cally brighter than the low-mass ones by a factor of ~2 and ~3 
for Class II and Class III, respectively. 



' There is significant but not complete overlap between CTTS and 
Class II stars and between WTTS and Class Ills. In particular, among 
the 57 CTTSs, 2 are classified as a Class I, 3 1 as Class II, 14 as Class III 
and 10 have no nIR classification; among the 110 WTTSs, one is a 
Class I, 29 are Class II, 70 are Class in, and 10 are have no nIR classifi- 
cation. 



5. Results 

For the stellar samples introduced in §|4]we first present the re- 
sults of our non-parametric analysis of variability amplitudes as 
a function of timescale and rotational phase. Next, as a tool to 
constrain the origin of variability, we investigate the incidence 
of flaring activity on the same samples. 

5.1. Variability amplitudes 

As described in §[3] we defined the representative variability 
amplitude, VA(At), of a given source sample and for a given 
timescale At, as the 50% width (diff'erence between 75% and 
25% quantiles) of the distiibution of all the F(t)/F(t + At) that 
can be computed from the X-ray light curves of stars in that 
sample. We now discuss how this quantity depends on At and 
on stellar characteristics by comparing diff'erent source samples. 
Next, we repeat the same comparative analysis as a function of 
stellar rotational phase (instead of timescale) by normalizing the 
time-lag At by the stellar rotation period, Pr„i, i.e. by considering 
VAiA(l>) = VAiAtlPro,). 

5.1 .1 . Timescale 

The left-hand column of Fig.|4]shows the run of VA vs. timescale 
in the interval 0.4-11 days for CTTS, WTTSs, Class II, and 
Class III stars. The three rows refer to the three X-ray energy 
bands: full, soft, and hard (from top to bottom). The contribu- 
tion of Poisson statistics to the VAs, as determined from Monte 
Carlo simulations, is also shown (both median and Icr disper- 
sion). Note that this latter contribution is not to be subtracted 
linearly from the observed VAs. In the hypothesis the distribu- 
tions of \0gF\IF2 are Gaussian (which is not far from being 
true, cf. Fig.|2]i, the Poisson contribution should be subtracted in 
quadrature. They thus contribute little to most of the measured 
VAs. The left-hand panels in Fig.|4]lead to the following conclu- 
sions: 

1. All samples are significantly variable on all timescales. 

2. Variability (VAs) increases strongly with timescale, up to at 
least 1 1 days, i.e. the longest timescale probed by the COUP 
observation. However, for some samples and energy bands 
(in particular CTTSs in the soft band), VAs may saturate at 
long timescales. 

3. CTTS and Class II stars appear to be systematically more 
variable than WTTS and Class III stars. 

4. The difference in VA between CTTS/ClassII and 
WTTS/Class III stars seems to increase toward long 
timescales. This is especially true for CTTSs and WTTSs in 
the soft band. 

We also studied the dependence of variability on stellar mass. 
Figure [5] shows again, repeated from Fig. [4] the run of VA vs. 
timescale for the two richest samples. Class II and Class III stars, 
along with the trends obtained for the two mass-segregated sub- 
samples of stars in the same classes. We see that, while the trends 
for the mass-segregated subsamples are qualitatively the same 
as for the full samples, low-mass stars of both classes are more 
variable than their higher mass counterparts. 

5.1.2. Rotational phase 

The left-hand column of Figure|4] shows the run of VA vs. the 
rotational-phase diff'erence A(p, ranging from 0.05 to A.QP,ot- 
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Fig. 4. Variability amplitude, VA, as a function of timescale (left column) and rotational phase difference (right column), for CTTS, 
WTTS, Class II, and Class III stars (see legends within panels), and for three different X-ray energy bands: full (0.8-8.0 keV, top 
row), soft (0.5-1.5 keV, middle), and hard (1.5-8.0keV, bottom). As described in §[3] the VA for each sample and x-axis value (At 
or A0) is defined as the 50% width of the distribution of \0gF1IF2 (cf. Fig.|2]), i.e. of the difference between the 75% and 25% 
quantiles. Formal error bars reflect photon statistics and are estimated thorough Monte Carlo simulations. Thick lines refer to the 
VAs derived from the data, while thin lines with error bars and matching color indicate the contribution of Poisson noise (also 
gstimated with Monte Carlo simulations). 
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Fig. 5. Variability amplitude, VA, as a function of timescale for 
Class II and Class III, for the full samples (same as in Fig.|4|l and 
for two mass-segregated subsamples (see caption). 



The purpose of these plots is to investigate the presence of ro- 
tational modulation, as evidenced by peaks at A(p - 0.5 Piot and 
A(p = 1.5 Plot- Other possible peaks at higher A<p values (e.g. 
2.5) cannot be investigated with our data since in order to have 
significant number of stars in each A(p bin we had to choose pro- 
gressively wider bins with increasing A(p. We also note that, in 
contrast to the case of the VA{At) plotted on the left, the num- 
ber of stars entering in the definition of the VA(A<p), i.e. those 
that define the distributions of log F((p)/F{<f) + A(p), generally de- 
creases with For example, given the fixed length of our X- 
ray observation, the variability at higher values can only be 
probed in stars with shorter rotational periods. From these plots 
we draw the following conclusions: 

1 . Several samples show, in one or more spectral bands, signs 
of rotational modulation. In particular small but statistically 
significant peaks at A(p - 1.5 and/or 0.5Prot are observed in 
all three bands for WTTS and Class III stars. 

2. Less clear signs of rotational modulations are generally seen 
for CTTS and Class II stars, except in the soft band, and in 
particular for CTTSs, for which a clear modulation signal is 
revealed by pronounced peaks at both 0.5 and 1.5Piot. The 
peaks are significant at the 3.8 and 2.3cr level, when com- 
pared to the mean of the two adjacent point^and suggest ap- 
parently larger amplitudes with respect to WTTSs/Class III 
stars. 



5.2. Flare frequencies and lengths 

To quantify the flaring activity for the same samples of CTTSs, 
WTTSs, Class II, and Class III stars discussed above, we applied 
the direct-detection method discussed in §3.2 Figure|6] shows 
the rate of flares as a function of minimum flare counts and the 



* For both VA(At) and VA(A<f)) the number of usable flux-ratios per 
star varies along the x-axis. 

' The significance of the peak at 1 .5 Ppoi would be higher if computed 
with respect to the minimum at 0.75 Piot- 



cumulative distribution of measured durations for our four main 
samples. We conclude that the flares from the four samples have 
very similar properties, both in terms of rate of occurrence and 
in their duration. 

We note, however, that this observation alone does not al- 
low us to conclude that the contribution of flares to the variabil- 
ity amplitudes, VAs, is the same for all samples. Indeed, VAs 
are a measure of relative variability while our flare rates indi- 
cate that the absolute intensity of flares is the same. However, 
since we have seen in §|4]that CTTSs (and, to a lesser extent. 
Class II stars) are on average fainter than WTTSs and Class III 
stars, the same absolute variability (i.e. flare rate) may translate 
into a more pronounced relative variability (i.e. VA) in CTTSs 
than in WTTSs, as is indeed observed. It remains to be seen, 
however, if this effect is strong enough to justify quantitatively 
the observed differences in VAs among our samples. 

In the following section we resort to Monte Carlo simula- 
tions to constrain the contribution of flaring activity to the vari- 
ability amplitudes as a function of timescale. 

6. Discussion 

We now use the results obtained in §|5]to discuss the origin of 
variability in the X-ray band. We consider two separate mecha- 
nisms: magnetic flares and rotational modulation. First, we try to 
understand to what extent the observed variability and the differ- 
ences we observe among subsamples can be explained in terms 
of flares. We then focus on rotational modulation and, finally, 
discuss the physical implications for the considered samples. 

6.1. Tlie contribution of flares 

The effect of flares on the fractional variability amplitudes de- 
pends, in addition to flare intensity and frequency, also on their 
durations and on the quiescent or characteristic count-rate on 
top of which they are observed. For example, for a flare of a 
given intensity (emitted energy or number of detected counts), a 
longer duration or a higher characteristic count-rates will result 
in lower relative variability amplitudes. Moreover, the trend of 
flare-induced VAs with timescale will depend on flare lengths: 
since flares are expected to occur stochastically with no par- 
ticular timescale, rising trends are only possible on timescales 
shorter than a few decay times, t. 

We have seen that our four main samples share similar fre- 
quencies and durations for detectable flares (§5.2i. However, 



since the samples have mean fluxes (or count-rates) that differ 
by up to a factor of ~ 2, (§|4]i, we cannot straightforwardly con- 
clude that the effect of flares on the VAs is the same. To assess 
the effect of flares on VAs realistically (even if not rigorously), 
we ran Monte Carlo simulations of light curves that attempt to 
reproduce the following observables: /) the time-averaged source 
count-rates, //) flare frequencies, and ///) flare durations. These 
simulations are fully described in Appendix [A] and lead to the 
conclusion that flares are likely not the dominant source of vari- 
ability as measured by the VAs and cannot explain the differ- 
ences between our samples. This is particularly true at long 
timescales. We indeed show that for realistic flares durations the 
run of VAs with timescale remains rather flat. 



6.2. Rotational modulation 

We have seen that the run of VA with phase lag shows signa- 
tures of rotational modulation, both for WTTS/Class III stars 
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Fig. 6. [left]: Rate of occurrence of flares as a function of minimum number of flare counts for our four main samples. The legend 
indicates the color of the lines representing each sample as well as the respective numbers of detected flares and sources, [right]: 
Cumulative distribution of flare durations for the same four samples. 



and, most strikingly, for CTTSs in the soft X-ray band. To inves- 
tigate this scenario in greater detail we computed the expected 
VA signal for a purely sinusoidal modulation with a given rela- 
tive amplitude A (Flux - 1 .0 -H A sin (p). 

For this simple experiment we drew rotational periods ran- 
domly from either a uniform distribution (between one and ten 
days) or from the actual distribution of rotational periods in 
our samples, with almost identical results. Rotational phases at 
the beginning of the observation were also chosen randomly. 
Figure]?] [left] shows the VAs vs. phase-lag trends expected for 
relative amplitudes ranging from 5% to 30%. For better compari- 
son with observations we then added the "modulated" VA signal 
to a constant VA level, meant to reproduce the contribution of 



all non-modulated variability sources such as flaring (see §6.1 



as well as Poisson noise. Because the modulated and constant 
contributions have to be added in quadrature, the absolute value 
of the latter smaller contribution has an important effect at 
values where the modulated VAs are low, but little effect on the 
peaks at Aip-O.S and 1.5 P™,. 

Figure]?] [right] compares the resulting simulated VAs with 
the observed trends for WTTS and CTTS in the low-energy X- 
ray band (lower and upper gray lines). Modulation amplitudes 
and constant contributions were chosen to reproduce the ob- 
served VAs at A0=O.O6 (the leftmost points) and A(p-Q.5 (i.e. 
where the first maxima are observed). For both WTTSs and 
CTTSs, very good agreement with the observed trends can be 
obtained for < 0.6? d, i.e. up to the datapoints following 
the peaks at 0.5 d, by adding in quadrature a constant contri- 
butions of ~0. 16dex and ~0.18dex and for modulation ampli- 
tudes of 18% and ~29%, respectively. For higher values of Acp 
the agreement is not as good: for CTTSs, even though the maxi- 
mum predicted at A0=1.5 Proi is indeed observed, the minimum 
predicted at A<^=1.0 Prot is too deep and the one at 2.0 Pni is not 
observed. Similar considerations apply to WTTSs. The quanti- 
tative disagreement may indicate that the actual variability is not 
sinusoidal, as assumed, and/or, maybe more importantly, that the 
X-ray modulation periods are not the same as the rotational pe- 



riods (see below), and/or that the structures responsible for the 
modulation (e.g. coronal active regions or absorbing circumstel- 
lar structures) evolve on timescales of a few rotational periods 
so that the modulation signal is lost for A0 > 1-2 fro/- Given the 
above considerations, we conclude that the clearer modulation 
signal observed for CTTSs in the soft X-ray band indeed sug- 
gests stronger modulation than for WTTSs. 

A similar analysis for the hard X-ray band, for which only 
WTTS show clear indications of rotational modulation, results 
in a best estimate for the modulation amplitude of WTTS of 24% 
and a non-modulated contribution to the VAs of ~0.22dex. For 
CTTS, instead, our simple model is not able to account for the 
observed trend, indicating that rotational modulation plays a rel- 
atively less important role with respect to other sources of vari- 
ability such as flares (which cannot account for the increase of 
VAs with timescale, however) and the evolution of emitting or 
absorbing structures (i.e. active regions and circumstellar mat- 
ter). 

Our detection of rotational modulation in the soft band is 



fairly consistent with the conclusions of Flaccomio et al. ( 2005 i, 
who found signatures of rotational modulation in 23 out of the 
233 COUP sources in their parent sample, with X-ray amplitudes 
in the full energy band ranging from 20% to ?0%. No evidence 
was found of a correlation between accretion and the detection of 
X-ray modulation and/or modulation amplitudes, apparently at 
odds with our findings that suggest larger modulation amplitudes 
for CTTSs, at least in the soft X-ray band. The small sample 
size of| Flaccomio et aL] ( 2005| l, differences in the samples and 
the fact that only the full X-ray band was considered might very 
well explain the potential discrepancy. 



Our samples significantly overlap with that of jFlaccomio 
|et al.| ( 2005) , who adopted for their analysis a lower count-rate 
threshold with respect to ours and excluded stars with short ro- 
tational periods, Pk,, <2d. In our sample of CTTSs (WTTSs), 
for example, 1? (51) out of the 19 (?2) stars in our "full band" 
sample are in common between the two studies. Among these. 
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Fig. 7. [left]: VA vs. rotational phase for simulated lightcurves. The simulations reproduce sinusoidally varying stars, observed with 
the COUP observation and analyzed in the same way as the real data. Four curves show results for modulation amplitudes from 
5% to 30% , as labeled, [right]: Comparison of the trends observed in the soft band for WTTS and CTTS (green and red lines, 
respectively) with two simulated VA curves (labeled gray lines) obtained by summing in quadrature the VAs for pure sinusoidal 
modulation (with amplitudes of 18% and 29%) with constant non-modulated VAs (0.16 and 0.18 dex, see labels). Modulation 
amplitudes and constant contributions are chosen to reproduce the observed VAs at A0=O.O6 (the leftmost points) and A0=O.5 (the 
first maxima). 



|Flaccomio et al.'(2005) claimed the detection of rotational mod- 
ulation in three CTTSs and ten WTTSs. 

Although |Flaccomio et al. ( |2005[ ) did not detect a difiference 
between CTTS and WTTS with regard to the incidence of X- 
ray rotational modulation or modulation amplitudes, a hint was 
found at the 1% significance level that CTTSs (and Class II stars) 
have X-ray modulation periods that are preferentially one half 
of the optically determined periods. The same hint is present in 
the subsample in common with the present study, since all three 



CTTS for which a periodicity was detected by Flaccomio et al. 
( 2005 1 were found to have Px iay ~ 0.5 Popt, while the same is 
true tor only two out of ten WTTSs. If, as it seems to be the 
case for CTTSs, X-ray periods are shorter than the optically de- 
termined period in a significant fraction of cases, our interpre- 
tation of the peaks in the runs of VA vs. phase lag (these latter 
computed using Popt) will be affected. In the cases for which 
Px-iay ~ 0.5 Popt we would indeed expect VA peaks at 0.25 Popt, 
0.75 Popt, etc. Although these peaks would not be detectable with 
the A(p binning we adopted in our analysis (e.g. Fig. [4]), they 
might fill in the minima between the maxima expected for the 
Px-ray ~ Popt case (e.g. between the 0.5 Popt and 1.5 Popt). As 
already mentioned, this might explain why the observed peaks 
and minima, e.g., in the VAs for CTTSs in the soft X-ray band, 
are not as pronounced as predicted by our simulations of sinu- 
soidally modulated lightcurves (Fig.|7]i. 

6.3. Implications 

Building on the above discussion, we now speculate on the phys- 
ical mechanisms responsible for the observed X-ray variabil- 
ity, with particular reference to the difference between CTTSs 



and WTTSs, i.e. by our definition stars undergoing "strong" and 
"no" accretion, respectively. We focus on the following results: 
flares are a certain source of variability in all our subsam- 
ples; ;/) most of the variability, however, cannot be explained 
by flares of the observed duration, especially at long timescales; 
///) the higher variability of CTTS/Class 11 stars with respect to 
WTTS/Class III stars cannot be explained by different flaring ac- 
tivities; iv) flux modulation due to stellar rotation is a good can- 
didate to explain at least part of the non-flaring variability; v) in 
the soft X-ray band, in particular, CTTSs show significant sig- 
natures of rotational modulation, suggesting higher modulation 
amplitudes than in WTTSs; vi) no indication of rotational mod- 
ulation is observed for CTTSs in the harder X-ray band. 



6.3.1 . The coronae of CTTSs 

The enhanced rotational modulation observed for CTTSs in the 
soft X-ray band with respect to the hard band may indicate that 
the azimuthal distribution of observable structures emitting soft 
X-rays is more inhomogeneous than that of those responsible for 
the hard X-rays. Two scenarios are qualitatively consistent with 
this conclusion: ;) modulation of the X-ray emission from accre- 
tion spot(s) or from an accretion-fed corona (as proposed by, e.g. 
,Brickhouse et al.|2010 1, both expected to be predominantly soft, 
and ii) modulation of the line-of-sight extinction toward a sig- 
nificant fraction of the coronal structures, again affecting mostly 
soft X-rays. The first accretion hypothesis builds on the obser- 
vations of X-ray emission from accretion spots, including re- 
cent observations of its rotational modulation (e.g. Argiroffi et al. 
|201 1 2012 1. The second absorption hypothesis could be related 
to the periodic obscuration of large parts of the stellar surface 
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by dusty structures in the inner disk, or by the gaseous accretion 
streams in co-rotation with the star Obscuration of the photo- 
spheric emission by disk warps has been extensively observed 
for AA Tau (Bouvier et al.||1999l [2 003 ; Me nard et al.||2003| 
[Bouvier et al.|2007a, ,Grosso et al.|2007j an d has been inferred 
to be common in CTTSs by Alencar et al. ( 2010[ l. Absorption 
by accretion streams was predicted to be effective in the X-ray 
regime by |Gregory et al. ( 2007) . The same accretion and absorp 



tion hypotheses were recently also proposed by Flaccomio et al. 
(|2010l) to explain the correlation between the X-ray and optical 



variability in the NGC 2264 star-forming region, found only for 



CTTSs and in the soft X-ray band. Flaccomio et al. (2010) fa- 
vored the absorption/ AA Tau-like scenario mainly because the 
contribution of the accretion shock emission to both the optical 
and the soft X-ray band was estimated to be too low to explain 
the large observed variability amplitudes. 

The same argument applies here to the X-ray emission from 
accretion spots: in all cases this emission has been revealed by 
the observation of emission lines from cool and dense plasma, 
detected by the Chandra and XMM-Newton high-resolution 
spectrographs. Even for the most extreme case known, TW Hya, 
the accretion shock emission is too soft and too faint to con- 
tribute a large part of the broad-band X-ray flux detected by 



imaging instruments, such as the Chandra ACIS (Dupree et al. 
|20T2). 



Moreover, the presence of such a significant soft contribu- 
tion would be evident in the observed CCD spectra and would 
require, for CTTSs only, a cool thermal component in the X-ray 
spectral fits. The AA Tau-like hypothesis, i.e. the heavy absorp- 
tion of part of the emission from the gaseous accretion streams, 
might instead find support from the time-average CCD spectra if: 
i) two thermal components of similar temperaure, but subject to 
very different absorptions, were required to fit the spectra or, //) 
the absorptions derived from the fits to the X-ray spectra were 
significantly higher than those derived from the optical colors 
and spectral types. The former scenario would indicate the pres- 
ence of two spatially distinct plasma components of comparable 
emission measure, only one of which intercepts the additional 
absorber (disk warp or accretion stream) along the line of sight, 
while the emission of the other is only subject to interstellar ab- 
sorption. The latter scenario might instead result if the circum- 
stellar structures absorb the great majority of the emission. 

Unfortunately, none of the above predictions is borne out 
by the CCD spectra of the 11 (19) CTTSs and the 62 (72) 
WTTSs entering in our VA analysis in the soft (full) X-ray 
bancj^ However, we argue that while for the accretion scenario 
we should indeed expect a detectable soft emission component, 
if it were responsible for the observed rotational modulation, the 
non-detection of an absorption "anomaly" in CTTS is actually 
consistent with the absorption/ AA Tau-like scenario. We reach 
this conclusion through Monte Carlo simulations that aimed to 
reproduce the observation of X-ray spectra composed of two in- 
trinsically similar thermal components absorbed by different gas 



* As for the presence of a significant accretion-spot contribution, we 
compared the temperatures of the cool thermal components needed to 
fit the CCD spectra of CTTS and WTTS, as reported by |Getman et al.| 
( |200 5 1. We found that CTTSs host, if anything, systematically hotter 
plasma than WTTS, at odds with the accretion hypotesis. As for the 
possible evidence of AA Tau-like behavior within our low-energy sam- 
ple, we found no evidence of double absorptions, both as indicated by 
[Getman et al.| (2005l, who flagged spectral fits with evidence of this 
situation, and upon new examination of the spectral fits. We also found 
no suggestion that CTTS have a peculiar Nu-Ay relation with respect to 
WTTS. 



columns. These simulations, described in detail in Appendix[B] 
show that the detection of the heavily absorbed component is 
extremely difficult, if not impossible, using current-quality data. 
We therefore conclude that the CCD spectra disfavor the accre- 
tion scenario as an explanation of the observed X-ray variability, 
but are compatible with the absorption/ AA Tau-like scenario. 

6.3.2. Time-averaged activity levels of CTTS and WTTS 

The absorption scenario is appealing also because it may pro- 
vide an explanation not only for the X-ray variability, which is 
the main subjet of the present study, but also for the previously 
recognized differences in the time-averaged X-ray activity lev- 
els of CTTSs (Lx or Lx/Lhoi) when compared to WTTSs: their 
on-average lower and more scattered activity levels at any given 
mass, or Lb,,!. 

Obscuration of large parts of the coronal plasma by the dense 
accreting columns or disk warps may indeed explain the obser- 
vation that CTTSs have, on average, 2-3 times lower activity lev- 
els than WTTSs, even assuming that their coronae are similar. As 
shown in Appendix |B] the obscured plasma cannot be detected 
and accounted for with current-quality data. Moreover, for suf- 
ficiently high absorbing column densities, Nfj > 3 x 10^^ cm"^, 
and sufficiently extended absorbing structures, intercepting >l/2 
of the coronal plasma, the retrieved X-ray luminosities will in- 
deed be less than one-half of the full coronal luminosity, with no 
discernible effect on the Nh estimated from the spectral fits. The 
observed magnitude of our VAs, up to ~0.5 dex for timescales 
of ~10 days and maybe larger for longer timescales, if attributed 
to periodic obscurations of parts of a corona from circumstellar 
structures, imply that the obscured fraction is significant. 

The estimate for the modulation amplitudes obtained in § |6.2| 
for the soft X-ray band, ~30%, if attributed to absorption implies 
a decrease in average flux on the same order Therefore, to ex- 
plain a reduction of a factor of ~2 in the observed emission of 
CTTSs, on average, we would need to assume that an additional 
fraction, ~20%, of the emission is permanently obscured. This 
seems reasonable. 

We now consider the difference in the scatter of the observed 
activity levels, making the reasonable assumption that the time 
variability of the stars in our samples is uncorrected. Regardless 
of its physical origin, time variability on long timescales will 
therefore contribute to the observed scatter in activity levels at 
any given mass or Lhoi- Our finding that CTTSs are more vari- 
able with respect to WTTSs may therefore explain the observed 
differences at least qualitatively. 

We now aim to assess the contribution of variability to this 
scatter more quantitatively. Indeed, our measure of variability, 
the VAs, would translate directly into a contribution to the 50% 
scatte:[^of activity levels (strictly speaking averaged over 30 ks) 
if they could be measured on the longest timescales (with the 
reasonable assumption that they level off at some timescale). 
Since we only measured variability amplitudes for timescales up 
to ~10 days, we consider our maximum VA for a given sample, 
yA(lOd), as a lower limit to the long-term variability, VA(oo). 

yA(lOd) is 0.47 dex for CTTS and 0.34 dex for WTTSs (cf. 
the upper-left panel of Fig.|4|i. We may compare these numbers 
with the scatter in the characteristic Lx vs. Lhoi relation found 
by Preibi sch etaL] ( |200y ) for larger samples of CTTS and WTTS 
in the COUP dataset, where the definition of WTTS and CTTS 



is the same as adopted here. Preibisch et al.| (j2005Jl, however. 



' The difference between the 25% and the 75% quantiles of the resid- 
uals. 
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Table 2. Parameters of best-fit linear regressions between Lx 
(three different sets of values, see text) and Lboi or mass (for 
M<0.8 Mq), and resulting 50% scatter. 





X 


a [erg 






b 


50% scatter [dex] 


Y 


CTTS 


WTTS 


CTTS 


WTTS 


CTTS 


WTTS 




^bol 


29.4 


30.1 


0.5 


1.1 


1.18 


0.65 




^bol 


29.7 


30.3 


0.4 


1.0 


1.18 


0.63 




^bnl 


29.5 


30.1 


0.4 


1.0 


1.06 


0.58 


^X.iOk.^ 


Mass 


30.0 


30.9 


1.4 


2.2 


1.07 


0.57 


Lxmv 


Mass 


30.2 


31.0 


1.4 


2.0 


1.04 


0.53 


^Xxlmr 


Mass 


30.0 


30.8 


1.7 


2.0 


0.99 


0.43 



only report the Icr dispersions (0.72 dex and 0.52 dex for CTTS 
and WTTS, respectively). To compute 50% dispersions for di- 
rect comparisons with our VAs, we repeated the Lx vs. Li,,,! re- 
gression analysis with the COUP sources using three sets of X- 
ray luminosities: the characteristic values, Lx^dmr, the time aver- 
ages Lx^av (over the 850 ks COUP exposure), and those relative 
to the individual 30 ks segments we employed for our variability 
analysis, Lz,30A.!- These latter values should allow the most direct 
comparison between VAs and scatter, and were derived by scal- 
ing the time-averaged Lx to the ratio between the count-rate in 
each 30 ks segment and the average count-rate. Throughout the 
regression analysis we used maximum likelihood (x^) linear fits 
that, by definition, yield minimum scatter 

Our regression analysis for the Lx.dmr vs. Lhoi relation of 
WTTS and CTTS yields very similar results to those presented 
by |Preibisch et al. (2005 i (c.f their Fig. 17). Best-fit parameters 
and 50% dispersions for these regressions and for those com- 
puted using Lx,av and Lx,3oiis are listed in Table|2] The same ta- 
ble also shows the results of regression analyses of the trends of 
Lx with mass for stars with masses below 0.8 Mq, i.e. for fully 
radiative stars, as shown in Fig.|8](for L^.c/mj)- Since the disper- 
sions around the maximum likelihood relations are systemati- 
cally lower (by 0.07-0.15 dex) than for the relation with Lhoi, we 
took these lower values as the measure of the true scatter of ac- 
tivity levels. 

When comparing these observed scatters in Lx with our 
VAs, we must acknowledge that the VAs are computed from 
smaller (and biased) subsamples with respect to the scatters. 
Fig. [8] shows with filled and empty symbols the stars belonging 
and not belonging to the count-rate limited subsample used for 
VAs (> Icount/Ks, or, roughly, logLx > 29.5 ergs"'). In the 
following we assume that the VAs of X-ray-faint stars, which 
also have statistically lower Lhoi and mass, are similar to those 
of their brighter counterparts. 

All dispersions listed in the last two columns of Table|2]for 
CTTSs and WTTSs, and in particular those for Lx^ioks which 
should be more directly comparable to the VAs, are higher than 
the observed yAs(lOd): 0.34 dex for WTTS and 0.47 dex for 
CTTSs. In performing this comparison, however, we should con- 
sider two facts that might explain this differenece: first, as men- 
tioned above, our VAs are computed for subsamples that, com- 
pared to the full samples, are biased toward high masses and 
bolometric luminosities. An extrapolation of our results on the 
dependence of VAs on mass (Fig.|5]l suggests that low-mass 
stars not included in our variability study may have higher VAs, 
and indeed Fig. [8] may suggest an increase of the scatter at low 
masses. Therefore, our VAs might be underestimated if applied 
to the whole sample down to the lowest masses and Lhoi- Second, 
and maybe most importantly, part of the observed scatter in the 
correlations between Lx and Lh„i or mass will be due to sta- 
tistical uncertainties in both quantities, e.g. introduced by the 
correction of the observed X-ray flux for the individual stellar 



absorptions, and/or by the derivation of stellar masses through 
placement of the star in the theoretical HR diagram and compar- 
ison with PMS evolutionary models, and/or by the inclusion in 
the samples of unresolved binaries. Following the discussion in 
|Preibisch et al.| ( |2005| l, we estimate that the 50% scatter intro- 
duced by the above uncertainties on the Lx vs. mass relation is 
~0.43 de?!*"! This contribution is to be subtracted (in quadrature) 
from the observed scatters before comparison with yA(lOd). 

Given these considerations, we thus argue that for WTTSs, 
a significant part of the scatter in Lx at a given mass might in- 
deed be justified by variability on the timescale of ~10 day^ 
The scatter for CTTSs is, however, significantly larger than the 
corresponding VA, even considering uncertainties, and cannot 
be explained in this way. Variability on longer timescales might, 
however, be responsible for most of the scatter in this case. A 
forthcoming study of longer-term variability will aim to address 
this point. In the absorption (or AA Tau-like) scenario, which we 
invoked to explain the evidence of significant rotational modu- 
lation in the soft X-ray band, this long-term variability might 
be naturally attributed to the year-over-year evolution of circum- 
stellar structures (disk warps and accretion streams), and the en- 
suing variability in the amplitude of the shading of coronal struc- 
trures (both modulated and unmodulated). This kind of evolu- 
tion is, indeed, well documented for AA Tau, and may also be 
inferred from the observed evolution of the structure of the mag- 
netic field in V2129 Oph ( ,Donati et al.,2010; . 



7. Summary and conclusions 

We have performed a study of X-ray variability among the young 
stars in the ONC, as observed with the 850 ks long COUP obser- 
vation. By investigating variability for different stellar samples 
as a function of timescale and rotational phase (for stars with 
known rotational period), we have found a number of novel re- 
sults: 

- CTTSs are significantly more variable than WTTSs at all 
timescales and in all X-ray spectral bands. 

- For all samples variability amplitudes increase with in- 
creasing timescale at least up to ~ 10 days, i.e. the longest we 
probed. 

- Signs of rotational modulation are visible, most clearly for 
CTTSs and in the soft 500-1500 keV band. 

- Among low-mass Class II and Class III stars, variability de- 
creases with increasing stellar mass. 

We speculate that the difference in variability between 
CTTSs and WTTSs may be explained assuming that the X-ray 
emission of CTTS is affected by time-variable absorption caused 
by circumstellar structures, such as warps in the inner disk and/or 
accretion streams, in co-rotation, or close to co-rotation, with 
the star. This suggestion is appealing because, assuming that the 
coronae of CTTSs and WTTSs are similar, it may also explain 
why CTTSs have lower and more scattered X-ray emission lev- 
els than WTTSs. 



Acknowledgements. We thank Jeremy Drake for stimulating discussions and 
careful reading of the manuscript. 



Derived summing in quadrature, 0.2 dex (uncertainties on masses), 
0.15 dex (Lx), and 0.15 dex (unresolved binaries), and converting this 
Icr scatter to a 50% scatter assuming that the distributions are Gaussian. 
" This is also confirmed by the scatter measured in the Lx.dmr vs. 
mass relation, i.e. using the X-ray luminosity averaged over lOd, and 
after excluding flares, which is on the same order as the estimated mea- 
surement uncertainties. 
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Fig. 8. Characteristic Lx vs. stellar mass for WTTSs and CTTSs (left and right panels, respectively). The samples are selected in 
the same way as for our variability study except that in this case we did not impose an X-ray flux (or count-rate) limit. Filled circles 
refer to the count-rate-limited sample used for the variability study. Linear fits are computed from all points (both filled and empty 
symbols) within the x-axis range spanned by the gray lines that represent them. 
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Appendix A: Monte Carlo flare simulations 



As discussed qualitatively in the main text (§ 6.1 1, the effect of 



flares on the fractional variability amplitudes, the VAs, depends 
on the intensity, frequency, and duration of flares, as well as on 
the quiescent or characteristic count-rate on top of which they 
are observed. In this appendix we employ Monte Carlo simu- 
lations to assess, quantitatively and realistically, the contribution 
of flaring activity to the VAs. We are particularly interested in de- 
termining whether flares can explain the observed rising trends 
of VAs with timescale and the differences among our four main 
samples (CTTSs, WTTSs, Class II, and Class III stars) . 

Although these samples share similar frequencies and dura- 
tions (§5.2 1 for detectable flares, our flare detection algorithm 
can only identify flares that individually contribute a statistically 
significant signal to the lightcurves. The VAs may, however, be 
affected also by smaller flares, either individually or by their su- 
perpositions. To assess the effect of flares on VAs realistically, 
we ran simulations of light curves that attempt to reproduce the 
following observables: /) the time-averaged source count-rates, 
//) flare frequencies, and ///) flare durations. More specifically, 
following the results of Caramazza et al. (2007 1, we assumed 
that the X-ray emission of our sources is produced by a super- 
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position of flares with impulsive rise phases and exponential de- 
cays, with a given fixed exponent t (at least initially). This last 
assumption is an approximation, since the observed flares clearly 
show a range of decay times, and was eventually relaxed. 

Flares are assumed to occur at an energy-dependent rate, 
dRfi(E), with a power law dependence: dRfi{C)/dC oc C^", 
where we substituted the energy with the number of detected 
flare counts, C. jCaramazza et aL] ( |2007 ) found a 2.2 for low- 
mass ONC members. We adopted the same value here, which is 
well-compatible with the bright tail of the flare frequency distri- 
bution in Fig. [6] (the turn-off' at log C < 3 may be explained by 



the inefficiency of the flare-detection process). Caramazza et al. 



( |2007 ) was able to reasonably reproduce the mean characteris- 
tics of observed light-curves, as well as the flare frequency dis- 
tribution, by assuming that the whole stellar X-ray emission is 
caused by flares, with a source-dependent minimum flare ampli- 
tude, Cmin, and total flare frequency, Rfi= J dRfi(C). Note that, 
if a < -2, a given set of a, C„,i„ and R ji, directly determines the 
mean source count-rate, as well as the characteristic one. 

Here, we followed a slightly simplified approach with re- 
spect to i Caramazza et aL]p007) : instead of trying to reproduce 
the characteristics of each individual light curve, i.e. by choosing 
individual Cm,„ and Rfi values to fit the mean and characteris- 
tic fluxes, we represented all the stars in each of our samples 
with a single set of Cmm and R fi values and thus tried to repro- 
duce the median count-rate of the stars in the sample. Moreover, 
in the last two of the four simulation sets we discuss below, we 
also allowed for the presence of a quiescent non-flaring emission 
component to better reproduce the observations, and, in the last 
set, for a range of flare decay timescales, t. 

All our simulations were run with a - -2.2 and o- = -2.4 but 
we only discuss the case of a = -2.2 since the conclusions are 
the same for the two cases. We chose C„„„ and Rfi so that /) the 
expected frequency of simulated flares with log C > 3 matches 
the observed value in Fig.|6j ~0.4 flares/Msec for all samples, 
and ii) the simulated lightcurves match the median count-rate 
of our samples. In particular we ran simulations to reproduce 
median count-rates of 3, 4, and 6cnts/Ksec (cf. Fig.|3]l. 

For each simulation set, i.e. for each set of parameters, a, 
Cmin, Rfi, and T, we generated 1000 simulated lightcurves, cov- 
ering the same temporal observing windows as in the COUP ob- 
servation (see e.g. Fig[T]i, and applied the exact same analysis as 
previously described for the real COUP sources in our samples, 
i.e. the VA analysis as a function of timescale and flare detection. 

For our first set of simulations we set r^lOks, i.e. similar to 



that of most observed flares. The top row of Fig. A.l shows the 
results. From left to right: the run of VAs vs. timescale, the flare 
frequency Rfi vs. flare counts, and the distribution of flare dura- 
tions. In all cases the black lines show the actual observed func- 
tions/distributions for our four samples, repeated from Figs. [2] 
and |6l The three thin red lines instead show the results of the 
simulations, one for each of the median count-rates, 3, 4, and 
6 cnts/Ksec, that encompass the observed one^^ We conclude 
that our models, while reproducing the observed flare durations 
reasonably well, produce light curves that are much more vari- 
able and with many more faint flares than observed. A possible 
solution to reduce both VAs and the frequencies of detected faint 
flares is to increase r. The second row in Fig. |A.l| shows that for 
T=100ksec we indeed obtain a good match with the observed 
VAs and flare frequencies. For VAs, we even qualitatively match 
the observed rising trends with timescale. However, the simula- 



tions do not match the measured flare durations, a discrepancy 
that is clearly supported by simple inspection of the observed 
and simulated lightcurves. 

In our third simulation set we therefore decided to revert t to 
10 ksec and, to lower VAs and flare frequencies, to substitute the 
flux emitted by small flares below a given counts threshold, with 
a quiescent emission. The third row in Fig. A. 1 shows the re- 
sult of setting the minium flare-counts threshold to 300. We note 
that /) flare frequencies are now reasonably close to the obser- 
vations, especially considering that our sharp cut of flare counts 
at 300 is obviously unrealistic; //) VAs are significantly smaller 
than observed, especially for long timescales, which leaves room 
for other sources of variability such as rotational modulation (of 
which we found evidence in § |5.1.2| i; Hi), flare durations are on 
average similar to the observations, even though with a signifi- 
cantly smaller scatter, as would be expected with our simplifying 
assumption of a single r value. To correct this single disagree- 
ment, in our last set of simulations (fourth row in Fig. A. 1 1 we 
allowed t to span a range of values, log-normally distributed 
with a mean=10ks and cr=l dex. We are thus able to reproduce 
the observed distribution and range of measured flare lengths, 
while the VAs and flare rates remain very similar to those of our 
previous simulations set. 

From the last two sets of simulations, the most successful at 
reproducing the observed flare statistics, we conclude that flares 
are likely not the dominant source of variability, as measured by 
the VAs. This is particularly true at long timescales. In fact, in 
our last simulation set with a reasonable mix of flare decay-times 
that well reproduce the measured flare durations, the run of VAs 
with timescale remains rather flat, which is not compatible with 
the observed sharp positive trends. 



Appendix B: Spectral Monte Carlo simulation 

Our preferred interpretation of the higher and more modulated 
X-ray variability of CTTSs with respect to WTTSs is the inter- 
vention of time-variable absorption of part of the X-ray emitting 
plasma. In this scenario the observable X-ray spectra of CTTSs 
would be the sum of two intrinsically similar components, un- 
dergoing very different attenuations: while the emission of one 
plasma component would be absorbed by interstellar material 
only, the emission from the other spatially distinct component 
would also intercept much thicker circumstellar material. This 
hypothesis seems to be contradicted by the time-averaged X-ray 
spectra obtained by the COUP project, showing neither signa- 
tures of such composed spectra, nor evidence of significantly 



higher extinctions for CTTSs (§ 6.3 i 



Note that in the plot of VA vs. timescale, higher VAs, a measure of 
relative variability, are associated with lower count-rates. 



To verify the significance of this lack of evidence, we have 
investigated through Monte Carlo simulations whether the two 
above signatures may or may not actually be detectable with 
data of the same quality as those in our hands. More specifically, 
we simulated the expected X-ray spectra in the two-absorptions 
scenario and verified whether the composite nature of the result- 
ing spectra could be detected with the customary X-ray spectral 
analysis. 

The simulated Chandra ACIS spectra were produced using 
XSPEC vl2.7 (Arnaud| |l996| l, adopting an isothermal model 
(APEC) with kT=1.26keV, quite typical of PMS stars, and as- 
suming zero interstellar extinction. The normalization was cho- 
sen so to obtain 1800 and 18000 counts (in two separate simu- 
lation sets). The lower value is somewhat higher than the mean 
and the median number of counts obtained by COUP for CTTS 
(as defined in this work): ~450 and ~1500, respectively. The 
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Fig.A.l. Results of our VA and flare analysis on simulated lightcurves consisting of a superposition of flaring and continuous 
emission (see text). From left to right we show measured VAs, flare frequencies, and flare durations. In all cases the three red 
lines refer to simulation sets that reproduce a given median count-rate: 3, 4, and 6 cnts/Ksec, spanning the median count rates of 
our four main stellar samples. Observed quantities for these four samples are drawn in black, reproduced from Figs. |2] and |6] The 
four rows refer to four different simulation sets that differ for the assumed flare decay-times and, in last two cases, for a lower 
cutoff in flare intensities (in counts), with a corresponding addition of a continuous emission level. From top to bottom: T=10ks (no 
cutoff/continuum), T=100ks (no cutoff/continuum), r=10ks (>300 counts plus continuum), and T=10ks±l dex (> 300 counts plus 
continuum) . 
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higher value is higher than the number of counts collected for all 
CTTSs, but one. 

For each simulation set, we simulated "IT, Iwo-Nh" com- 
posed spectra by varying /) fabs, the fraction of the coronal 
plasma that is absorbed (from 0. 1 to 1 .0, in 0. 1 steps), and //) Nh, 
the column density of the absorbing material, (from 10^' cm"^ to 
3.2 X lO^"* cm"^, in 11 logarithmic steps). For each point in this 
grid, we simulated 100 ACIS spectra, regrouped them so to have 
SNR~5 per spectral bin, and fit them, within XSPEC, with an 
usual two-temperature model with a single absorption ("2T, one- 
Nh")- We judge the ability to detect the two absorptions from the 
goodness of these latter fits, and specifically from the mean value 
of the null probability, (i.e. the probability of obtaining a higher 

value given the model), n.p.: high n.p. values indicate that a 
model with a single Nh is a good representation of the composite 
spectrum and that the absorbed component cannot be detected. 
To reject the single Nh hypothesis, we would require low values 
of n.p., e. g. <5 % (a conservative threshold). 

Figure B.l is a contour plot showing n.p. derived from the 
simulations with the brighter template spectrum (18000 counts 
for the unabsorbed case) in the fabs-Nn plane. We see that for 
much of the plane n.p. is larger than 5%, indicating that even 
with high statistics we are unable to detect the existence of 
two absorption components for most fab.rNH combinations. The 
same plot for the case of the fainter template spectrum (1800 
counts) shows that n.p. >10% for all fabs-Nn pairs. The mid- 
dle panel in Fig. |B.l| shows the mean unabsorbed flux obtained 
from the fits with 2T, one-Nn spectral models, at each position 
in the grid and relative to the flux for the Nh=0 case. To obtain 
a reduction of the flux by a factor of 2-3, i.e. on the order of 
the mean Lx difference between WTTS and CTTSs, Nh must be 
greater than ~ 3x10^^ cm"^ and fabs-O-5-Q.l . The contour lines 
relative to a decrease in flux by a factor of 2 and 3 are also re- 
ported in white in the n.p. contour plot (upper panel). Note that 
there is ample room in the fabs-Nn plane for having a decrease 
in retrieved flux on this order, and at the same time being unable 



to detect the second absorption. The bottom panel of Fig. B.l 
shows the mean Nh obtained from the spectral fits, showing that 
in the same region of the parameter space discussed above, not 
only is it impossible to detect the second Nh, but the value of 
the best fit Nh obtained from a single-Nn model is very low, ba- 
sically that of the unabsorbed component. We also note that in 
the same region of the plane, the values of the two temperatures 
obtained from the spectral fits with 2T, one-Nn models, are in 
line with observed values (not shown). 

We did not try to fit the spectra with more than two thermal 
components or to allow plasma abundances to vary in the spec- 
tral fits. Therefore the combinations of fibs-Nn values for which 
a single-A^H model is not adequate to fit the two-Nn spectra may 
be even more limited than found above. Moreover, we neglected 
systematic uncertainties, e.g. on the models and on the calibra- 
tion of the spectral responses, which, if considered, might make 
the rejection of a one-Nn model even harder 

We conclude that the scenario in which a large part of the 
coronal emission is obscured by stable or time-varying gaseous 
structures is compatible with the observed spectra, showing nei- 
ther evidence of double absorptions nor substantially higher best 
fit Nh values. Moreover, assuming that WTTS and CTTS have 
similar coronae and that the difference in observed flux is due to 
partial obscuration of the corona, we must require that the ob- 
scuring material has Nh greater than a few 10^^ cm"^. 
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Fig. B.l. Results of the simulations of IT, two-Nn spectra 
(-18000 counts for Nh = 0) fitted with 2T, one-A^^ models. 
The black contours in the upper panel show the average null 
probability as a function of the fraction of obscured emission 
measure, fibs, and the column density of the obscuring mate- 
rial, Nh- The middle and lower panels show contour plots of 
the mean unabsorbed flux (in the 0.5-8.0 kV band, in units of 
1.3x10"'-' s^'cm"^), and of the mean Nh (units 10^^ cm"^), both 
obtained from the spectral fits. The two white contours in the upj 
per panel are repeated from the middle panel and show the loci 
corresponding to a 50% and 70% reduction in measured flux. 



